Learning and inference in a nonequilibrium Ising model with hidden nodes 
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We study inference and reconstruction of couplings in a partially observed kinetic Ising model. 
With hidden spins, calculating the likelihood of a sequence of observed spin configurations requires 
performing a trace over the configurations of the hidden ones. This, as we show, can be represented 
as a path integral. Using this representation, we demonstrate that systematic approximate inference 
and learning rules can be derived using dynamical mean-field theory. Although naive mean-field 
theory leads to an unstable learning rule, taking into account Gaussian corrections allows learning 
the couplings involving hidden nodes. It also improves learning of the couplings between the observed 
nodes compared to when hidden nodes are ignored. 
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I. INTRODUCTION 

Within the statistical mechanics community, a signifi- 
cant body of recent work deals with parameter learning 
in statistical models. This work typically focuses on the 
canonical setting of the inverse Ising problem: given a 
set of configurations from an equilibrium [l| , or nonequi- 
librium [2[ Ising model, how can we find the interactions 
between the spins? These inverse problems can be diffi- 
cult to solve using brute force approaches 3 and thus for 
practical purposes, the results of such theoretical works 
offering approximate learning and inference methods are 
of great importance. The recent interest in using statis- 
tical models to study high-throughput data from neural 
P-Q, genetic Q and financial [3] networks further em- 
phasizes the necessity of developing efficient approximate 
inference methods. 

Despite its practical relevance, however, little has been 
done on learning and inference in models with hidden 
spins using statistical physics. Learning the parameters 
of a statistical model in the presence of hidden nodes and 
inferring the state of hidden variables lie at the heart of 
difficult problems in statistics and machine learning. It 
is well known that for statistical modeling of data, hid- 
den variables most often improve model quality [l^ . 
Furthermore, when used for network reconstruction and 
graph identification, ignoring hidden nodes may cause the 
model to incorrectly identify /remove couplings between 
observed nodes when such couplings do not exist / actually 
exist. For these reasons, it is quite important to study 
models with hidden nodes, both from a theoretical per- 
spective and for applications to real data where usually 
only a fraction of the system is observeable [1, [E, 0| • 

In this paper, we study the reconstruction of the cou- 
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plings in a nonequilibrium Ising model with hidden spins 
and the inference of the state of these spins. Theoreti- 
cally, the couplings can be reconstructed by calculating 
and maximizing the likelihood of the data which involves 
performing a trace over the configurations of the hidden 
spins at all times: an impossible task even for moderately 
sized data sets and networks. This is, in fact, a general 
problem, not only for a nonequilibrium Ising model, but 
for any kinetic model with hidden variables. It is, how- 
ever, natural to use a path integral approach for doing 
this trace, and this is what we do here. We show that 
efficient approximate learning and inference can be de- 
veloped by evaluating the log-likelihood in this way and 
employing mean-field type approximations. 

In what follows, after formulating the problem, we first 
evaluate the likelihood of the data using the saddle point 
approximation to the path integral. We find that the 
result allows inference of the state of hidden variables if 
the couplings are known but cannot be used for learning 
the couplings. We then perform a Gaussian correction 
to this saddle point approximation to derive Thouless- 
Anderson-Palmer (TAP) like equations [llj for our prob- 
lem. Our numerical results show that successul learning 
of the couplings, even those that involve hidden nodes, 
can be achieved using the corrected learning rules. We 
find that with enough data and using these equations 
even hidden-to-hidden connections can be reconstructed. 

Although in this paper our focus is on the paradig- 
matic case of the Ising model, the approach taken here is 
quite general and can be adapted to many other popular 
statistical models. 



II. FORMULATION OF THE PROBLEM 

Consider a network of binary spin variables subject 
to synchronous update dynamics. Suppose some of these 
nodes are observed for T time steps, while the other nodes 
are hidden. We denote the state of the observed spins at 
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time t by s{t) = {si{t)} and the hidden ones by cr = 
Wa{t)}- For clarity, we always use the subscripts . . . 
for the observed and a,b, . . . for the hidden ones. We 
consider the following transition probability 

p[{s,fT}it + l)\{s,a}{t)]^ (la) 

cxp[^^ s,{t + l)g,{t) + <^a{t + l)ga{t)]Z{tr^ 

Z{t) = []2cosh[5,(i)] 2cosh[g„(0] (lb) 



where gi{t) and ga{t) arc the fields acting on observed 
spin i and hidden spin a at time t, respectively 



9^{t) ^ 
9a{t) 



E 



^ JySj(t) + J,bOb{t) (2a) 

Ja,,S, (t)+^^Jaba6(t), (2b) 

Jij, Jai, Jia and Jaf, are the observed-to-observed, 
observed-to-hidden, hidden-to-observed and hidden-to- 
hidden connections. These couplings need not be sym- 
metric and, thus, the network may never reach equilib- 
rium. Although, here we consider the case of zero exter- 
nal field, our derivations follow also for nonzero fields. 

Given this model, we would like to answer the follow- 
ing questions: What can we say about the state of hidden 
spins from the observations? How can we learn the vari- 
ous couplings in the network? 

Answering these questions requires finding the likeli- 
hood of the data. Optimizing the likelihood with respect 
to Js yields the maximum likelihood value of the cou- 
plings. The posterior over the state of hidden nodes given 
the observed nodes is 



p[{^(t)}Lil{«(i)}Li] = 



P[{sit)}I=^] 



which also requires calculating the likelihood. 

The likelihood of the observed spin configurations un- 
der the model Eq. [1] is 

p[{sit)}t,] = Tr^Y[p[{s, CT}{t + ms, crm (3) 



The log-likelihood of the data is recovered at V' — ^ 0. 
Furthermore, Clip] acts as the generating functional for 
the posterior over the hidden spins: ma{t), the mean 
magnetization of the hidden spin a at time t given the 
data, can be found using 



ma{t) 



lim pLa{ip,t), 



/^a('0,O 



(6) 



Higher-order derivatives of C with respect to '0 yield 
higher-order correlation functions of the hidden spins. 

To develop mean-field approximations, we first per- 
form the trace in Eq. [5] using the standard approach for 
dealing with a nonlinear action in the local fields [H, [13 '■ 
we write the term inside the log in Eq. [5] as a path in- 
tegral over gi{t) and ga{t) enforcing their definitions in 
Eq.[2]by inserting (5[g,(t)-X]j Jv^ji't) -J2b Jib<^b{t)] and 
^[9a{t)-Y.j Ja]Sj{t)-J2b Jab<yb{t)] in the integral. Using 
the integral representation of the delta functions, C can 
be written as 



(7a) 

MtWit)] (7b) 



C[ip] = log j D^cxpi^] 
<I> = logTrcr cxp[Q + A 



where Q = {gi, gi, ga, ga}- A, as well as gi and ga, come 
from the integral representation of the S{-) functions used 
to enforce the definitions in Eq. [5] 



(8) 



+ J2^9ait)[ga{t) - Y JajSj{t) - ^ Jab(Tb{t)\ 
a,t j b 

The crucial thing about this rewriting of C is that the 
resulting action $ is now linear in a and we can thus 
easily perform the trace over a. The ability to perform 
the trace here comes at the cost of the high-dimensional 
integral over Q in Eq. I7al The integral, however, can be 
approximated using the saddle point approximation and 
the corrections to the saddle point. 



and involves a trace over {(t(1), . . . ,a{T)} which is dif- 
ficult to do. 

To perform this trace analytically, we consider the fol- 
lowing functional 

m ^ logTr. n t}(< + l)\{s, ajit)]. 

t 

(4) 

which, using the definition of the transition probability 
from Eq. [TJ can be written as 



£[ip] = logTrc^exp 



Q[s,(t] 



a.t 



(5) 



i,t 



log2cosh[(7i(t)] 



log2cosh[gQ(t)] 



III. SADDLE-POINT APPROXIMATION 

The saddle point values of ga,9i,9a and gi can be de- 
rived by putting the derivatives of $ with respect to these 
variables to zero. For '0 = 0, the saddle point equations 
read 

igiit) = tanh[.gi(t)] - s,{t + 1) (9a) 
iga{t) = tanh[ga(t)] - ma{t + 1) (9b) 
g^{t) = -^^r^jl*) + Eb J^bmb{t) (9c) 
9a{i) = Yj J'^jSjit) + J'^brribit) (9d) 

in which ma(t) satisfy the self-consistent equations 
ma{t) = (10) 
tanh [ga{t - 1) - iY,.9j{t)J]a - » Eh 5fc(*)'^ba] 
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Let us take a moment to describe the physical meaning 
of these equations. The saddle point equations for and 
5a, Eq- l9cland[9dl are just the mean values of the local 
fields at the mean magnetizations ma{t) of the hidden 
variables (see Eq. [2]). Examining Eq.[9aland l9bl one can 
identify —igi{t) and ~iga{t) as errors back propagated in 
time. The self-consistent equations for ma, Eg. I10[ take 
the form of the usual naive mean-field equations except 
that back propagated errors are taken into account. 

Given a sequence of configurations from the observed 
spins, solving Eq. [TU] yields the mean magnetization of 
the hidden ones given this data. The derivatives of $ 
with respect to the Js, evaluated at the saddle point 
solutions, give the maximum likelihood learning rules for 
the couplings within this saddle point approximation. 



IV. GAUSSIAN CORRECTIONS 

The saddle-point equations and learning rules can be 
improved by Gaussian corrections, that is by evaluating 
the contributions to the path integral from the Gaussian 
fluctuations around the saddle [3 Ell • This will lead to 
what can be called as TAP equations for our problem. In 
what follows, by doing this on the Legendre transform of 



a,t 

and using the equation of state, 

dm 



(11) 



we correct the naive mean-field equations for our problem 
with the hidden spins. 

Denoting the saddle point value of C hy Cg, and per- 
forming a Legendre transform we have (see Appendix IX)) 

ToM -5I[*^(* + l)5'(*)-log2cosh(5,(t))] (13) 

i,t 

+ ^[/ia(t + l)9a{t) - l0g2 COsh(5a(0)] + ^ S[fla{t)] 
a.t a^t 

where S'[/ia(t)] is the entropy of a free spin with magne- 
tization iJ-ait), and gi and ga are as in Eq. [5] with ab(t) 
replaced by ^h(t). Using Eq. [T2] and putting -0 = 0, we 
recover Eq. [TOl as we should. 

The Gaussian corrections can be performed as follows; 
for details see Appendix[Bj We first perform the Gaussian 
integral around the saddle point of C. This modifies Cq 
to 



Ci=Co-^logdct[d^C] 



(14) 



where d^C is the Hessian, containing the second deriva- 
tives of C with respect to g and g, calculated at their 
saddle point values. Keeping terms quadratic in J we 



then derive the modified equations for rua (t) by first cal- 
culating dCi/dipait)- Expressing Ci in terms of the new 
niait), again keeping terms quadratic in J, and finally 
performing the Legendre transform, we get 

FiH = EoH - ^ - tanh2g,(t)]4[i _ ^2^^^^ 

i.a.t 

- J E["^a(i + 1) - tanh^ .9a(0].4\[l - <it)]- (15) 



a,b.t 



where now, and ga are as in Eq. [2]but with instead 
of CTa. Eq. [15] together with the equation of state at — > 
0, 



dniait) 



0, 



(16) 



yields the new self-consistent equation for ma{t). The 
couplings can then be found by optimizing Fi with re- 
spect to the couplings. 

In sum, the corrected TAP-learning will be the follow- 
ing EM like iterative algorithm: at each step, solve the 
equations for rria derived from Eq. [16] given the current 
value of the couplings, and then adjust the couplings pro- 
portional to dTi/dJ. 



V. NUMERICAL RESULTS 

We evaluated the saddle point and TAP equations on 
data generated from a network of ~ 100 spins with 
the dynamics in Eq. [1] and with the couplings drawn in- 
dependently from a Gaussian distribution with variance 
Ji/N. We first studied how much the saddle point Eq. 
[To] and TAP equations Eq. \W\ tell us about the state of 
hidden spins if we know all the couplings. Fig. [T] A and 
B show the distributions of ma{t) when 20% of the net- 
work is hidden. The solutions to the TAP equations are 
more pushed towards and —1 compared to the sad- 
dle point. Using a simple estimator aa{t) = sgn(mo(i)) 
to infer the state of hidden spin a and time Fig. [T]C 
shows that the two approximations perform equally well 
in inferring the state of hidden spins. The effect of TAP 
corrections becomes very important in reconstructing the 
couplings. We found that the saddle point learning was 
always unstable leading to divergent couplings. However, 
including the TAP corrections dramatically changed this. 
Fig. [2] shows an example of this when 10% of the spins 
were hidden. In this case, good reconstruction of all cou- 
plings could be achieved. 

Although for a proportionally small number of hid- 
den nodes we observed good reconstruction, when the 
hidden part was large, for a constant data length, we 
found that the algorithm starts to become unstable: the 
root mean squared error for the inferred hidden-to- hidden 
connections increases after a few iterations and the oth- 
ers follow. This instability, however, seems to arise from 
our attempt to learn the hiddcn-to-hidden connections. 
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FIG. 1. (Color online) Solutions for ma- (A) and (B) show 
histograms of ma{t) found from solving saddle point and TAP 
equations, Eqs. 1101 1161 given the correct couplings. Here we 
had Ji = 0.7, 20% hidden spins and a data length of T = 
10^. The dark gray bars show the histogram of the values for 
which aa{t) was actually +1 while the transparent ones show 
when it was actually —1. (C) shows the percent correct if 
we use aa{t) — sga{ma{t)) to infer the state of aa{t) vs. the 
percentage of the visible part. TAP and saddle point results 
were virtually indistinguishable. 
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FIG. 2. (Color online) Reconstructing the full network. Scat- 
ter plots showing the inferred couplings (blue dots) versus the 
correct ones in a network of A'' = 100 total spins, with 10 hid- 
den, Ji — 1 and data length of T = 10® . Red crosses show the 
initial values of the couplings used for learning. A, B, C, and 
D show observed-to-observed, observed-to-hidden, hidden-to- 
observed and hidden-to-hidden couplings, respectively. 



VI. DISCUSSIONS 



For an architecture where there is no hiddcn-to-hidden 
couplings, that is a dynamical semi-restricted Boltzmann 
machine [16|, or when hidden-to- hidden couplings exist 
in the original network but are ignored during learning 
(i.e. put to zero), the algorithm allows recovering the 
other connections even for very large hidden sizes. This 
is shown in Fig. [S] 

In the examples above, we have used the right number 
of hidden nodes when learning the couplings. In practice, 
however, the number of hidden nodes may be unknown. 
We thus wondered if the number of hidden nodes can be 
estimated from the data. In Fig. |31 we plot the objective 
function Fi/iV for the data used in Fig. [2l where N is 
the total number of visible nodes, that is, the sum of the 
observed nodes, and an unknown number M of hidden 
nodes. As can be seen in Fig. 21 varying M, at the correct 
number of hidden units a peak in Ti/N is seen and also 
the reconstruction error of the visible-to-visible connec- 
tivity reaches a minimum. This was the case, even when 
we constrained the hidden-to-hidden connections to zero 
during learning. However, the peak in the normalized 
cost function was sharper, and the minmum achieved for 
the reconstruction lower, when hidden-to-hidden connec- 
tions were also learned. 



In this paper we described an approach for learning 
and inference in a nonequilibrium Ising model when it is 
partially observed, that is, when the history of spin con- 
figurations are known only for some spins and not for the 
others. Treating the log-likelihood of the data as a path 
integral over the fields acting on hidden and observed 
spins in a non-equilibrium Ising model, we approximated 
it first using saddle point approximations, and then by 
considering TAP corrections to it. 

In usual statistical mechanics settings, TAP equations 
can be derived in several ways: using Plefka expan- 
sion 17 1, information- geometric arguments [isj or by 
calculatin g t he Gaussian fluctuations around the sad- 
dle point [Tsj . Here we took an approach similar to 
the latter case except that we did not use the Hubbard- 
Stratonovich transformation. It would be quite interest- 
ing to see what the other approaches yield as TAP correc- 
tions. Tyrcha and Hertz have also studied the prob- 
lem of hidden spins, using a different approach. Their 
equations agree with ours at the saddle point level but 
differ at the TAP level. 

In terms of performance, saddle point equations gave 
results that were similar to TAP for inferring the state 
of hidden spins, but were very unstable for learning 
the couplings. Using TAP equations, learning all the 
couplings, even the hidden-to-hidden ones, was possible 
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FIG. 3. (Color online) Root mean squared error of the 
couplings. (A)-(C) show the rms error for the observed- 
to-observed, hidden-to-observed and observed-to-hidden cou- 
plings using TAP learning when the original model did not 
have hidden-to-hidden connections; the other couplings were 
drawn independently from a Gaussian. (D)-(F) show the 
same but when the hidden-to-hidden connections were present 
in the original network but were ignored (i.e. put to zero) 
when learning the other couplings. In (A) and (D) the green 
dots show the inferred couplings by maximizing the likeli- 
hood with hidden nodes completely ignored. Here we used 
Ji = 1,A'^ = 100, T = 10^, and five different random realiza- 
tions of the connectivity. 



when the hidden part was proportionally small. By ig- 
noring hidden-to-hidden couplings, even when they were 
present in the network, we could achieve good perfor- 
mance in recovering all the other couplings, even with 
half the network hidden. This also improved recovering 
the observed-to-observed couplings compared to when 
hidden spins were ignored all together. Our numerical 
results also indicate that even the number of hidden units 
can be estimated by optimization of the objective func- 
tion. This is of particular practical revelance as often 
the true size of the network is unknown. In Fig. 21 both 
the objective function and the reconstruction quality in- 
crease as we add hidden nodes suggesting that even a 
wrong number of hidden nodes can improve the recon- 
struction of observed-to-observed couplings. This can be 
explained by the fact that the added hidden variables, 




• • Full 




FIG. 4. (Color online) Estimation of the number of hidden 
nodes. (A) the RMS error of the observed-to-observed cou- 
plings inferred using TAP learning for different values of as- 
sumed hidden units, M. The couplings are inferred either 
assuming no hidden-to-hidden connectivity (green x) or by 
trying to learn all couplings (blue circles). (B) the value of 
the objective function per inferred node for different numbers 
of assumed hidden units. The data used here is the same as 
that of Fig. [21 that is in reality there were 10 hidden nodes 
which is where the normalized objective function shows max- 
imum and the RMS error shows a minimum. 



even if wrong in number, to some degree can explain the 
correlations induced by the hidden spins in the real net- 
work. 

An important future direction, in our view, is to study 
how averaging the likelihood over some couplings, in par- 
ticular the hard to infer hidden-to-hidden ones, given a 
prior over them, will influence the reconstruction of the 
other couplings. There is a large body of work where 
the path integral approach is used to study the disorder- 
averaged behavior of disordered systems [H, [131 • The 
approach that we developed here can be easily combined 
with this work to average out some of the couplings. 
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Appendix A: Legendre transform 



Here we derive Eq. [T3]iii the main text. Starting from 
Eq. 7, wc have 



£[ip]=\og J Dgexp[<i>], (I-la) 
* = ^ s,it + l)g,{t) + i J2 Ut) [9^{t) - J2 J^r^At)] + * J2 [gait) - J2 Ja^Mt)] - log2 cosh[5,(t)] 

i,t i.t j a,t i i.t 

- ^I0g2c0sh[ga(t)] + ^ log 2 cosh [ga{t - 1) - i^Jbagtit) - i^Jiagiit) + ^a(t) 



9$ 



= tanh [ga{t - 1) - i'Ygj{t)Jja - i'Ygb{t)Jba + V'a(i)]; 



(I-lb) 
(I-lc) 



with the saddle point being at Eq. |9l with /.ia instead of point value of £o- Given Eq. Il-lbl we have 
rua- Applying the following identity to the last sum in 
Eq. intl 



log2cosh[.T] = 5'[tanh[a;]] + .Ttanh[x] 



S[x] 



1 + a; 



■log 



1 + a; 



1 - X 



log 



1-x 



(1-2) 
(1-3) 



and using Eq. II-lcl and the saddle point values for gi and 
(7a, we obtain 



*0 = Y[s^{t + l)g,{t) - Iog2c0sh[(7,;(i)]] 

+ + l)9a{t) - log2cosh[5a(t)]] 

+ 5^/.a(i)V'a(i)+E^[/^-Wl (1-4) 



At the saddle point we have £0 = 'i'o and using To 
Co - t '^a{t)^ia{t) we arrive at Eq. [I3l 



B 



tt' 



dg,{t)dg,{t') 



-%^tt'[l - tanh^[5i(t)]] 



r<tt 

'-'nh 



E, 



92$ 



dga{t)dgb{t') 

= -5ab5tt' \jjlit + 1) - i&nh^[ga{t)]] 

r 



F, 



tt' 



dga{t)dgb{t') 
92$ 



-iJibll - Aft(t)](5t_i,t' 



92$ 



i5ab5tt' -iJba[^ - MQ(i + l)]'5t',t+i 



and d'^^/dgi{t)dgj{t') — iiSy^tf leading to 
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Appendix B: Corrections to the saddle-point 

lilcelihood (11-2) 

as part of the Hessian for t and t' = t + 1, and the same 
structure gets repeated for other times. The first to the 
To include the Gaussian corrections in Eq. [TH]we first fourth rows/columns correspond to derivatives with re- 
calculate the Gaussian fluctuations around the saddle spect to gi, gi, ga, {ja all at t and the fourth to eighth 
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rows / columns at t' = t + 1. 

Denoting the matrix of the blocks on the diagonal part 
of the Hessian as a and the rest as /3, wc write 

log det(a + /3) = log det(a) + log det[/ + a"^/3] 

= logdet(a) + Trlog[/ + a-i/3] (II-3) 

« logdet(a) + Tr[a-^/3] + ^Tr{[a-^ P]^}} + ■■■ 

Assuming that the couplings are random and indepen- 
dent with mean of order and a standard deviation 
of Ji/^/N as typically assumed in mean- field models of 
spin glasses l2l|, logdet(a) will be of quadratic order in 
Ji. Since a is block diagonal, so will be and therefore 
Ti[a~^(3] = 0. Ignoring the second trace as contributing 
higher order terms, we thus approximate the determi- 
nant of the Hessian matrix by considering only the block 
diagonal terms 

det{a) = Y[det{A{t)B{t)+l} det{C{t)D{t)+l} (H-4) 



leading to the following estimate of the contribution of 
the Gaussian fluctuations 



dC^-lY.^ogdet{A{t)Bit) + l} 



■logdct{C{t)D{t) + 1} 



(H-5) 



^ E - tanh2[5.(i)]] E [1 - W]} 

td b 

4 E {[^'(i + 1) " tanh2(ffa W)] E -^aJl - (^im 



With these Gaussian fluctuations taken into account, the 
corrected value for the mean magnetization is 



ma(t) 
Ut) = 



_ dC 
dSC 



= M.W+'aW (H-6) 
(H-7) 



We now express Tq to the quadratic order in J in terms 
of nia instead of Ha- With some algebra at V we 
have 



ToH = J2^s,{t + l)g,{t) - log2cosh(.g,(i))] 

i,t 

+ Y^[mait + l)gait) - log 2 cosh(5,(t))] + J2 S[mait)] 

a,t a,t 

-^[si{t + 1) - tanh(.g,(t))] J„LW 

i,t.a 

+ Y,la{t + l)[tanh-iK(i + 1)] - gait)] 

a,t 

-Y,bMt + 1) - ta.nh[ga{t)]Jabh{t) 



.t,b 



Where we have abused the notation here, using ga and gi 
to indicate the fields calculated at nia (t) instead of fia (t) ■ 
Noting that la itself is quadratic in J, the terms that 
involve la above are all higher order, and can therefore 
be ignored. This, combined with the expression for SC 
in Eq. IH- 5 1 yield Eq. [T5]in the main text. 
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